{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "rbm_ansatz.ipynb",
      "version": "0.3.2",
      "views": {},
      "default_view": {},
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  },
  "cells": [
    {
      "metadata": {
        "id": "yJC-_S9dXh_U",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# RBM as a Quantum Wave Function Ansatz"
      ]
    },
    {
      "metadata": {
        "id": "wi81QOR46AVX",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "try:\n",
        "    import torch\n",
        "except:\n",
        "    !pip install torch torchvision\n",
        "    import torch\n",
        "import numpy as np\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "egWqWgGp6AVV",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Variational Monte Carlo\n",
        "With RBM as a state ansatz, we need to train its parameters in order to make it have maximum overlap with the ground state wave function of a hamiltonian.\n",
        "\n",
        "Notice that the cost function now is the energy rather than the log-likelihood with some known dataset:\n",
        "$$E_\\theta=\\frac{\\langle\\psi_\\theta|H|\\psi_\\theta\\rangle}{\\langle\\psi_\\theta|\\psi_\\theta\\rangle}$$\n",
        "Here, $\\theta$ represents parameters in a RBM.\n",
        "\n",
        "We use Variational Monte Carlo (VMC) method to get the energy expectation and gradients to optimize the ansatz.\n",
        "VMC states that the gradient of $E$ with respect to $\\theta$ is defined as \n",
        "$$g=\\langle E_{\\rm loc}\\Delta_{\\rm loc}\\rangle- \\langle E_{\\rm loc}\\rangle\\langle\\Delta_{\\rm loc}\\rangle$$\n",
        "where $\\langle\\cdot\\rangle$ is the ensemble average over Monte Carlo importance sampling $x\\sim|\\psi_\\theta(x)|^2$, $E_{\\rm loc}=\\frac{\\langle x|H|\\psi_\\theta\\rangle}{\\langle x|\\psi_\\theta\\rangle}$, $\\Delta_{\\rm loc}\\equiv\\frac{\\partial\\log\\psi_\\theta(x)}{\\partial \\theta}$."
      ]
    },
    {
      "metadata": {
        "id": "J3qdADxY6AVf",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### VMC measurements\n",
        "1. obtain some samples using `vim_sample` method.\n",
        "2. perform measurements on samples.\n",
        "\n",
        "\n",
        "VMC gives expectation values for $\\langle E_{\\rm loc}\\rangle$, $\\langle\\Delta_{\\rm loc}\\rangle$ and $\\langle E_{\\rm loc}\\Delta_{\\rm loc}\\rangle$.\n",
        "It requires a kernel with following methods\n",
        "* **propose_config**: propose a new configuration, given old configuration.\n",
        "* **prob**: give probability for a spin configuration $\\langle\\sigma|\\psi\\rangle$."
      ]
    },
    {
      "metadata": {
        "id": "YnySuxew6AVg",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "def vmc_sample(kernel, initial_config, num_bath, num_sample):\n",
        "    '''\n",
        "    obtain a set of samples.\n",
        "\n",
        "    Args:\n",
        "        kernel (object): defines how to sample, requiring the following methods:\n",
        "            * propose_config, propose a new configuration.\n",
        "            * prob, get the probability of specific distribution.\n",
        "        initial_config (1darray): initial configuration.\n",
        "        num_bath (int): number of updates to thermalize.\n",
        "        num_sample (int): number of samples.\n",
        "\n",
        "    Return:\n",
        "        list: a list of spin configurations.\n",
        "    '''\n",
        "    print_step = np.Inf  # steps between two print of accept rate, Inf to disable showing this information.\n",
        "\n",
        "    config = initial_config\n",
        "    prob = kernel.prob(config)\n",
        "\n",
        "    n_accepted = 0\n",
        "    sample_list = []\n",
        "    for i in range(num_bath + num_sample):\n",
        "        # generate new config and calculate probability ratio\n",
        "        config_proposed = kernel.propose_config(config)\n",
        "        prob_proposed = kernel.prob(config_proposed)\n",
        "\n",
        "        # accept/reject a move by metropolis algorithm (world's most famous single line algorithm)\n",
        "        if np.random.random() < prob_proposed / prob:\n",
        "            config = config_proposed\n",
        "            prob = prob_proposed\n",
        "            n_accepted += 1\n",
        "\n",
        "        # print statistics\n",
        "        if i % print_step == print_step - 1:\n",
        "            print('%-10s Accept rate: %.3f' %\n",
        "                  (i + 1, n_accepted * 1. / print_step))\n",
        "            n_accepted = 0\n",
        "\n",
        "        # add a sample\n",
        "        if i >= num_bath:\n",
        "            sample_list.append(config_proposed)\n",
        "    return sample_list\n",
        "\n",
        "\n",
        "def vmc_measure(local_measure, sample_list, measure_step, num_bin=50):\n",
        "    '''\n",
        "    perform measurements on samples\n",
        "\n",
        "    Args:\n",
        "        local_measure (func): local measurements function, input configuration, return local energy and local gradient.\n",
        "        sample_list (list): a list of spin configurations.\n",
        "        num_bin (int): number of bins in binning statistics.\n",
        "        meaure_step: number of samples skiped between two measurements + 1.\n",
        "\n",
        "    Returns:\n",
        "        tuple: expectation valued of energy, gradient, energy*gradient and error of energy.\n",
        "    '''\n",
        "    # measurements\n",
        "    energy_loc_list, grad_loc_list = [], []\n",
        "    for i, config in enumerate(sample_list):\n",
        "        if i % measure_step == 0:\n",
        "            # back-propagation is used to get gradients.\n",
        "            energy_loc, grad_loc = local_measure(config)\n",
        "            energy_loc_list.append(energy_loc)\n",
        "            grad_loc_list.append(grad_loc)\n",
        "\n",
        "    # binning statistics for energy\n",
        "    energy_loc_list = np.array(energy_loc_list)\n",
        "    energy, energy_precision = binning_statistics(energy_loc_list, num_bin=num_bin)\n",
        "\n",
        "    # get expectation values\n",
        "    energy_loc_list = torch.from_numpy(energy_loc_list)\n",
        "    if grad_loc_list[0][0].is_cuda: energy_loc_list = energy_loc_list.cuda()\n",
        "    grad_mean = []\n",
        "    energy_grad = []\n",
        "    for grad_loc in zip(*grad_loc_list):\n",
        "        grad_loc = torch.stack(grad_loc, 0)\n",
        "        grad_mean.append(grad_loc.mean(0))\n",
        "        energy_grad.append(\n",
        "            (energy_loc_list[(slice(None),) + (None,) * (grad_loc.dim() - 1)] * grad_loc).mean(0))\n",
        "    return energy.item(), grad_mean, energy_grad, energy_precision\n",
        "\n",
        "\n",
        "\n",
        "def binning_statistics(var_list, num_bin):\n",
        "    '''\n",
        "    binning statistics for variable list.\n",
        "    '''\n",
        "    num_sample = len(var_list)\n",
        "    if num_sample % num_bin != 0:\n",
        "        raise\n",
        "    size_bin = num_sample // num_bin\n",
        "\n",
        "    # mean, variance\n",
        "    mean = np.mean(var_list, axis=0)\n",
        "    variance = np.var(var_list, axis=0)\n",
        "\n",
        "    # binned variance and autocorrelation time.\n",
        "    variance_binned = np.var(\n",
        "        [np.mean(var_list[size_bin * i:size_bin * (i + 1)]) for i in range(num_bin)])\n",
        "    t_auto = 0.5 * size_bin * \\\n",
        "        np.abs(np.mean(variance_binned) / np.mean(variance))\n",
        "    stderr = np.sqrt(variance_binned / num_bin)\n",
        "    print('Binning Statistics: Energy = %.4f +- %.4f, Auto correlation Time = %.4f' %\n",
        "          (mean, stderr, t_auto))\n",
        "    return mean, stderr"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "nY9-FD6m6AVk",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### VMC Kernel definition\n",
        "When proposing a new configuration, it has 5% probability to flip all spin, can making VMC sample better in Heisenberg model, proposed configuration satisfies spin conservation.\n"
      ]
    },
    {
      "metadata": {
        "id": "gkG0kni56AVl",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "class VMCKernel(object):\n",
        "    '''\n",
        "    variational monte carlo kernel.\n",
        "\n",
        "    Attributes:\n",
        "        energy_loc (func): local energy <x|H|\\psi>/<x|\\psi>.\n",
        "        ansatz (Module): torch neural network.\n",
        "    '''\n",
        "    def __init__(self, energy_loc, ansatz):\n",
        "        self.ansatz = ansatz\n",
        "        self.energy_loc = energy_loc\n",
        "\n",
        "    def prob(self, config):\n",
        "        '''\n",
        "        probability of configuration.\n",
        "\n",
        "        Args:\n",
        "            config (1darray): the bit string as a configuration.\n",
        "\n",
        "        Returns:\n",
        "            number: probability |<config|psi>|^2.\n",
        "        '''\n",
        "        return abs(self.ansatz.psi(torch.from_numpy(config)).item())**2\n",
        "\n",
        "    def local_measure(self, config):\n",
        "        '''\n",
        "        get local quantities energy_loc, grad_loc.\n",
        "\n",
        "        Args:\n",
        "            config (1darray): the bit string as a configuration.\n",
        "\n",
        "        Returns:\n",
        "            number, list: local energy and local gradients for variables.\n",
        "        '''\n",
        "        psi_loc = self.ansatz.psi(torch.from_numpy(config))\n",
        "\n",
        "        # get gradients {d/dW}_{loc}\n",
        "        self.ansatz.zero_grad()\n",
        "        psi_loc.backward()\n",
        "        grad_loc = [p.grad.data/psi_loc.item() for p in self.ansatz.parameters()]\n",
        "\n",
        "        # E_{loc}\n",
        "        eloc = self.energy_loc(config, lambda x: self.ansatz.psi(torch.from_numpy(x)).data, psi_loc.data)\n",
        "        return eloc, grad_loc\n",
        "\n",
        "    @staticmethod\n",
        "    def propose_config(old_config, prob_flip=0.05):\n",
        "        '''\n",
        "        flip two positions as suggested spin flips.\n",
        "\n",
        "        Args:\n",
        "            old_config (1darray): spin configuration, which is a [-1,1] string.\n",
        "            prob_flip (float): the probability to flip all spins, to make VMC more statble in Heisenberg model.\n",
        "\n",
        "        Returns:\n",
        "            1darray: new spin configuration.\n",
        "        '''\n",
        "        # take ~ 5% probability to flip all spin, can making VMC sample better in Heisenberg model\n",
        "        if np.random.random() < prob_flip:\n",
        "            return -old_config\n",
        "\n",
        "        num_spin = len(old_config)\n",
        "        upmask = old_config == 1\n",
        "        flips = np.random.randint(0, num_spin // 2, 2)\n",
        "        iflip0 = np.where(upmask)[0][flips[0]]\n",
        "        iflip1 = np.where(~upmask)[0][flips[1]]\n",
        "\n",
        "        config = old_config.copy()\n",
        "        config[iflip0] = -1\n",
        "        config[iflip1] = 1\n",
        "        return config"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "VaEcgJmtVHdB",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Wave function ansatz\n",
        "Comparing with the model defined in last example, we make serveral adaptions here\n",
        "* using {-1, 1} valued nodes instead of {0, 1} valued nodes.\n",
        "* use probability function of visual nodes directly, $p(v)\\sim e^{x^T a}\\prod\\limits_i2\\cosh(W_ix+b_i)$."
      ]
    },
    {
      "metadata": {
        "id": "g_6uhS8x7fdH",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "class RBM(nn.Module):\n",
        "    '''\n",
        "    Restricted Boltzmann Machine\n",
        "\n",
        "    Args:\n",
        "        num_visible (int): number of visible nodes.\n",
        "        num_hidden (int): number of hidden nodes.\n",
        "\n",
        "    Attributes:\n",
        "        W (2darray): weights.\n",
        "        v_bias (1darray): bias for visible layer.\n",
        "        h_bias (1darray): bias for hidden layer.\n",
        "    '''\n",
        "\n",
        "    def __init__(self, num_visible, num_hidden):\n",
        "        super(RBM, self).__init__()\n",
        "        self.W = nn.Parameter(torch.randn(num_hidden, num_visible).double() * 1e-1)\n",
        "        self.v_bias = nn.Parameter(torch.zeros(num_visible).double())\n",
        "        self.h_bias = nn.Parameter(torch.randn(num_hidden).double() * 1e-1)\n",
        "        self.num_visible = num_visible\n",
        "        self.num_hidden = num_hidden\n",
        "        \n",
        "    def psi(self, v):\n",
        "        '''\n",
        "        probability for visible nodes, visible/hidden nodes here take value in {-1, 1}.\n",
        "\n",
        "        Args:\n",
        "            v (1darray): visible input.\n",
        "\n",
        "        Return:\n",
        "            float: the probability of v.\n",
        "        '''\n",
        "        v = v.double()\n",
        "        if self.W.is_cuda: v = v.cuda()\n",
        "        res = (2 * (self.W.mv(v) + self.h_bias).cosh()).prod() * (self.v_bias.dot(v)).exp()\n",
        "        return res"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7_O7XUu_6AVp",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Hamiltonian\n",
        "Here, we use the Heisenberg model with rotated basis\n",
        "$$H=\\sum\\limits_i J_z S_i^zS_{i+1}^z-\\frac{J_{xy}}{2}(S^+_{i}S^-_{i+1}+S^-_{i}S^+_{i+1})$$\n",
        "Comparing with traditional Heisenberg model, we performed transformation $S^{x,y}_i=-S^{x,y}_i,i\\in \\rm odd$, intended to make the ground state wave function always positive.\n",
        "\n",
        "Instead of defining the hamiltonian directlt, we define a function to calculate $E_{\\rm loc}$"
      ]
    },
    {
      "metadata": {
        "id": "7U5HEs2m6AVq",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "def heisenberg_loc(config, psi_func, psi_loc, J=1.):\n",
        "    '''\n",
        "    local energy for 1D Periodic Heisenberg chain.\n",
        "\n",
        "    Args:\n",
        "        config (1darray): bit string as spin configuration.\n",
        "        psi_func (func): wave function.\n",
        "        psi_loc (number): wave function projected on configuration <config|psi>.\n",
        "        J (float): coupling strengh. Here, J = Jxy = Jz.\n",
        "\n",
        "    Returns:\n",
        "        number: local energy.\n",
        "    '''\n",
        "    # get weights and flips after applying hamiltonian \\sum_i w_i|x_i> = H|x>\n",
        "    nsite = len(config)\n",
        "    wl, flips = [], []\n",
        "    # J*SzSz terms.\n",
        "    nn_par = np.roll(config, -1) * config\n",
        "    wl.append(J / 4. * (nn_par).sum(axis=-1).item())\n",
        "    flips.append(np.array([], dtype='int64'))\n",
        "\n",
        "    # J*SxSx and J*SySy terms.\n",
        "    mask = nn_par != 1\n",
        "    i = np.where(mask)[0]\n",
        "    j = (i + 1) % nsite\n",
        "    wl += [-J / 2.] * len(i)\n",
        "    flips.extend(zip(i, j))\n",
        "\n",
        "    # calculate local energy <\\psi|H|x>/<\\psi|x>\n",
        "    acc = 0.\n",
        "    for wi, flip in zip(wl, flips):\n",
        "        config_i = config.copy()\n",
        "        config_i[list(flip)] *= -1\n",
        "        eng_i = wi * (psi_func(config_i) / psi_loc)\n",
        "        acc += eng_i\n",
        "    return acc"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "rB2orpQedo1c",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## The training"
      ]
    },
    {
      "metadata": {
        "id": "KMzFhYgc6AVc",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "def train(model, learning_rate, use_cuda):\n",
        "    '''\n",
        "    train a model.\n",
        "\n",
        "    Args:\n",
        "        model (obj): a model that meet VMC model definition.\n",
        "        learning_rate (float): the learning rate for SGD.\n",
        "    '''\n",
        "    initial_config = np.array([-1, 1] * (model.ansatz.num_visible // 2))\n",
        "    if use_cuda:\n",
        "        model.ansatz.cuda()\n",
        "\n",
        "    while True:\n",
        "        # get expectation values for energy, gradient and their product,\n",
        "        # as well as the precision of energy.        \n",
        "        sample_list = vmc_sample(model, initial_config, num_bath=200, num_sample=8000)\n",
        "        energy, grad, energy_grad, precision = vmc_measure(model.local_measure, sample_list, num_spin)\n",
        "\n",
        "        # update variables using steepest gradient descent\n",
        "        g_list = [eg - energy * g for eg, g in zip(energy_grad, grad)]\n",
        "        for var, g in zip(model.ansatz.parameters(), g_list):\n",
        "            delta = learning_rate * g\n",
        "            var.data -= delta\n",
        "        yield energy, precision"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "e24315mq6AVv",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "### Start training\n",
        "We will see GPU is not better than CPU in this case. There are mainly two factors causing this\n",
        "* the scale of vectorized operation is not large enough\n",
        "* freqent data transfer between CPU and GPU can become a bottleneck."
      ]
    },
    {
      "metadata": {
        "id": "BEpsTZ5N6AVx",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          },
          "base_uri": "https://localhost:8080/",
          "height": 2249
        },
        "outputId": "59da2a02-b64c-4700-d010-9836ccca97f3",
        "executionInfo": {
          "status": "error",
          "timestamp": 1525697698874,
          "user_tz": -480,
          "elapsed": 55825,
          "user": {
            "displayName": "刘金国",
            "photoUrl": "//lh3.googleusercontent.com/-lDAT81T3HSE/AAAAAAAAAAI/AAAAAAAAAgw/eH3JEob7M1Y/s50-c-k-no/photo.jpg",
            "userId": "116824001998056121289"
          }
        }
      },
      "cell_type": "code",
      "source": [
        "%matplotlib inline\n",
        "import time\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# set random number seed\n",
        "use_cuda = False\n",
        "seed = 10086\n",
        "torch.manual_seed(seed)\n",
        "if use_cuda: torch.cuda.manual_seed_all(seed)\n",
        "np.random.seed(seed)\n",
        "max_iter = 100\n",
        "num_spin = 8\n",
        "num_hidden = 16\n",
        "E_exact = -3.65109341\n",
        "\n",
        "# visualize the loss history\n",
        "energy_list, precision_list = [], []\n",
        "def _update_curve(energy, precision):\n",
        "    energy_list.append(energy)\n",
        "    precision_list.append(precision)\n",
        "    if len(energy_list)%10 == 0:\n",
        "        plt.errorbar(np.arange(1, len(energy_list) + 1), energy_list, yerr=precision_list, capsize=3)\n",
        "        # dashed line for exact energy\n",
        "        plt.axhline(E_exact, ls='--')\n",
        "        plt.show()\n",
        "\n",
        "rbm = RBM(num_spin, num_hidden)\n",
        "model = VMCKernel(heisenberg_loc, ansatz=rbm)\n",
        "\n",
        "t0 = time.time()\n",
        "for i, (energy, precision) in enumerate(train(model, learning_rate = 0.1, use_cuda = use_cuda)):\n",
        "    t1 = time.time()\n",
        "    print('Step %d, dE/|E| = %.4f, elapse = %.4f' % (i, -(energy - E_exact)/E_exact, t1-t0))\n",
        "    _update_curve(energy, precision)\n",
        "    t0 = time.time()\n",
        "\n",
        "    # stop condition\n",
        "    if i >= max_iter:\n",
        "        break"
      ],
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Binning Statistics: Energy = -2.5752 +- 0.0455, Auto correlation Time = 0.4185\n",
            "Step 0, dE/|E| = 0.2947, elapse = 2.3380\n",
            "Binning Statistics: Energy = -2.6432 +- 0.0417, Auto correlation Time = 0.3898\n",
            "Step 1, dE/|E| = 0.2761, elapse = 2.4210\n",
            "Binning Statistics: Energy = -2.6183 +- 0.0472, Auto correlation Time = 0.6120\n",
            "Step 2, dE/|E| = 0.2829, elapse = 2.3998\n",
            "Binning Statistics: Energy = -2.6812 +- 0.0381, Auto correlation Time = 0.4503\n",
            "Step 3, dE/|E| = 0.2656, elapse = 2.4022\n",
            "Binning Statistics: Energy = -2.6423 +- 0.0393, Auto correlation Time = 0.5433\n",
            "Step 4, dE/|E| = 0.2763, elapse = 2.3924\n",
            "Binning Statistics: Energy = -2.7394 +- 0.0368, Auto correlation Time = 0.4653\n",
            "Step 5, dE/|E| = 0.2497, elapse = 2.4108\n",
            "Binning Statistics: Energy = -2.7198 +- 0.0313, Auto correlation Time = 0.3983\n",
            "Step 6, dE/|E| = 0.2551, elapse = 2.3942\n",
            "Binning Statistics: Energy = -2.8283 +- 0.0281, Auto correlation Time = 0.4056\n",
            "Step 7, dE/|E| = 0.2254, elapse = 2.4171\n",
            "Binning Statistics: Energy = -2.8373 +- 0.0328, Auto correlation Time = 0.5386\n",
            "Step 8, dE/|E| = 0.2229, elapse = 2.4173\n",
            "Binning Statistics: Energy = -2.8795 +- 0.0346, Auto correlation Time = 0.6192\n",
            "Step 9, dE/|E| = 0.2113, elapse = 2.4222\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VPWd//HXZHIjyeRGJgkJl9zg\nC+GmXAQEVNBi6w1KW2291P6Etrbd/f1aW3etq2tb2u3FXWu7XbcXbL2t1VrdWihWq7aoBFIBEQH5\nYoBwTzJgEgIhCUnm98ckJBMSIJOEmRzez8eDR2bO+Z45H7+dvvPN99xcfr8fERFxnqhwFyAiIgND\nAS8i4lAKeBERh1LAi4g4lAJeRMShosNdQDufr27Qn86TlpZAdXV9uMuIGOqPDuqLYOqPYH3pD6/X\n4+ppnUbw/Sg62h3uEiKK+qOD+iKY+iPYQPWHAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i\n4lAKeBERh1LAi4g41KAP+LsfKeHuR0rCXYaISMQZ9AEvIiLdU8CLiDiUAl5ExKEU8CIiDqWAFxFx\nKAW8iIhDDfqAb25p5XjDSezealr9g/6ZISIi/SZinugUqpPNrTQ0tfDDp98hPTmOGeOymFGcxYjM\nJFyuHh904ljt1wQ8+OVLw1yJiIRbSAFvjIkGHgUK2z7jG9bat7q0mdzWBuBFa+2yvhTak/hYN9Fu\nF5OKMthgq3ipdC8vle4lNyORmeOzmDEui4zUIQOxaxGRiBbqCP424Li1do4xZjzwG+CSLm1+CXwB\n2AT8jzEmwVrb7w9hdLlcxES7ueOacdy2YAzvlh2hdFsl7+48zPOrd/H86l0UDU9hVnEW08Zm4kmI\n7e8SREQiUqgB/xTw27bXPmBo55XGmCwgyVq7sW3RZ0LcT6/ERLuZNjaTaWMzqW84yXrro3RbJdv3\nVFO2v5anX/2A8fnpzByfxcVFXuJi9VxIEXGukALeWnsSONn29qvA012a5AEfGmMeA0YDz1lrHw6x\nxpAkxMdw2eQcLpucQ3VdI6XbKlm3rYLNO4+weecR4mLcTBmTwYzibMbnp+GOGvTHm0VEgpw14I0x\nS4GlXRY/YK192RjzFWAKcH2X9S4gH1gEnADWGmP+Yq3d2tN+0tISQnqyuNsdOJDq9Xp6bOP1ehhT\nkMFt141nX2UdqzfuZ/U7+1m7tZK1WytJSYpl7uRcLp8yHDMqLaSDs0u++woAj963oNfb9qdz6Y/z\nKdx1RMr/LhD+vog06o9gA9EfZw14a+1yYHnX5caYJQSCfVHbiL6zSmCrtfZIW9u3gPFAjwFfXR3a\n9HxLS+DUSJ+v7pzax0fB1dOGs2BqLrsOHmXd1kr+vr2SlWt2s3LNbryp8cwozmbW+CyGDU3sVR1u\nt+uc6xgozc2tuFzhrwMCX9hw19Hb78dAiYS+iCTqj2B96Y8z/WII9SyaAuBO4HJrbUPX9dba3cYY\njzEmHagBLiJw0DViuFwuCnNTKMxN4aYri9hWXk3ptgo27jjMypJyVpaUMzIriZnF2cwoziLNExfu\nkgHw+/0cb2jmcO0JDtc04Ks9ga+mgcM1J/DVNvBhXSNRLnjmtQ+YZjIpyE0m6gI8XVREQj/IupTA\ngdVVxpj2ZQuAu4DV1tq1wNeAlwA/8Gdr7bt9rLVb/XG+d7Q7ikmFQ5lUOJTGphbeKfNRurWSLbs/\n5Hd/LeO5v5ZhRqYyc3w204yXhPiYfqi8Z00nWzhc28DhtvD21ZwIvK85ga/2BCcaW7rdLmlIDO4o\nF62tfl55ex+vvL2P1KRYpppMphkvo4enEhWlsBe5UIR6kPVe4N5uVv2gU5tSYEaIdYVNXKybmcXZ\nzCzOpq6+ifXbq1i7rZLte2vYvreGp17ZweTCocwozmJy0VBiQjhu0Nrq58O6hsAIvG3kfWpEXnOC\n2uNN3W4XGxOFN2UIGcPjyUgdgjd1CN6UwOuMlHiGxEVz9yMl+P1+PvtRw/rtPt75wMdrG/bz2ob9\nJCfGMnWMl2nGy5iRqTqwLOJwg/5K1oHkSYhl3pThzJsynMM1Jyh9v5J1WyvZsMPHhh0+hsRFM9V4\nmVWchd/vJ3BsOTCNUnfiZGDkXdN1JH6CD4820tJ6+m0Volwu0pPjGDcqDW9qPBkpQ8hIjcebEghz\nT0LMOR0AdrlcTCrMYFJhBs0thu17q1m/3cfGHT7++s4B/vrOAZKGxDBljJdpY72MHZlGtFthL+I0\nLn+E3L/F56uLjELOwu/3s993nHVbK1i3rZLqukYAXC6IiY7CmzqEwzUNNJ7sfholOTEWb1toZ7SF\nePtIPC05rs+j6jPdqqCltZUde2tYbwO/oI62/aWQGB/NxaMDYV+cl95vYR+OA2nNLa0c8B2nvOIo\n5RV1rHmvgtjoKH78j7ND+murv+igYjD1R7A+HmTtcdSngO+DVr+fD/bVsG5bJW9sOogfGBLnDoy8\nU+IDwd02fdI+jRIXM7Ahc673omlt9fPB/rawt1XUHAuE/ZC4aC4qymD62EzG56eFHIp3P1KC2+3i\nB1+cFdL256KltZVDh+vZ3Rbm5Yfq2Fd1jOaW1tPapnniuHbWKOZOyiEm+vz/taJAC6b+CDZQAa8p\nmj6IcrkwI9MwI9N4b+cRoqJc/PDOWYPiJmdRUR21f+aq0ew6cJT1tor1toq1WytYu7WC+Fg3FxVl\nMNVkMrEgndgB/uV0Jq2tfio+rA+MzA/VUV5Rx97KOpqaO8LcHeViuDeJ/GEe8oYlk5ft4SfPbeZE\nUzPHG07y1Cs7WLVuD9fNymPOpGGalhLHU8D3E5fLRVSUa1CEe1dRLhdFw1MoGp7CTfOL2H2oLhD2\n26tYt62SddsqiYtxM6lwKNPGZjKpYOiA3ubB7/dTVX0iMDJvC/M9lXU0NnVMe0W5XORkJJI3zEN+\ndiDQh3sTT/uLIyrKRWJ8DPffPo2XSvfw+sYDPPGy5U9r93D97DwunZCtoBfHUsBLEJfLRUFOMgU5\nyXzqikL2Vh5jva3i7e0d/2Kjo5hYMJSpY71MLsxgSFzoXyO/38+R2gbKK+qCAv1EY3NHTcCwjETy\nsj2Bf8OSGZGZ1KvpruTEWG6aP5qrLxnJqnV7+Ns7B3nspe2sLCnnhtn5zJqQpbOKxHEU8A7Tn/eB\nd7lcjMr2MCrbw+LLCtjvO8767YFpnPYziaLdUUzIT2faWC8XFWWc8RoBv99PdV1jYL68U5gfOxF8\nIXRWegKTC4eeCvORWUnEx/bPVzU1KY6brxrDx2aMYtXaPax+9wC/XvU+K9eWc8PsPGYUK+jFORTw\nck5cLhcjMpMYkZnExy8r4IDvGOutj/W2ik1lh9lUdhh3lIvx+elMNd7A07VaYVPZYcoPtR0Erag7\ndeZOO29qPONGpZE3zENedjKjsjwkxA/81zLNE8ctC8bwsZkj+dPaPbzx7kGWr3yfFSV7WDg7j0vG\nZemiMBn0FPASklxvErneJBbOyefQkeOBs3G2V526W2e7n/5+86nXQ5PjmDrG2xHm2R6ShgzsVcFn\nk54cz21XGz42cyQrS/aw5r1D/HLFNlaUlLNwTj7TxmbqVg8yaOk0yX5yPk4LHAwqq+vZYH384c1d\nREW5uGbmKPKyA2e0JCdG/sNWfDUnWFFSTsl7FbT6/eR6E1k4O58pxtunoNdpgcHUH8F0HvwgoC9t\nh8H+C6+yup6Va8op2VqB3w/DvUksmpvPxaMzQjpTKhK+G5H0vN5I6I9IovPgRc6jrLQEllxXzLWX\n5rFizW7Wba3kZy+8x8isJBbNKWBy0dBBeUqsXFgU8CJnkJ2ewOevH8+1s/L445rdvP1+FT99fjN5\n2R4Wzc1nYoGCXiKXAl7kHORkJHLnwglcf+kxXlxTzvrtVTz83GYKcpJZNCef8fnpPQb9YJ+uksFL\nAS/SC7neJL68aAL7qo7xx7d2s2GHj4d+9y5FuSksnJtPcYiPfBxofr+fVr+fyKtMBpICXiQEIzKT\n+MriieypqOPFt3azqeww//HMJsYMT2HR3ALGjko7b7U0nmyh9lgj1XWN1BxrouZYY9u/JmrqOl63\n3+H0+dU7uWbmqD5dgSyDg86i6Uc6MyDYhdQfuw8d5cW3dp+6BmDsyFQWzS1gzIjUkKdomltaqe0S\n2NWnArsjwOs73dahKxfgSYwlNSmWQ0fqaWlppdUfuHXD4ssKmDNxWFgu6LqQvhvnQqdJDgL60ga7\nEPtj18Gj/OGtXWzZ9SEAxXlp7K86Tnyc+1TAt7b6OVrfFtx1nQO8kepO7+vquz7LPlhifDSpSXGk\neuJITYoNvG7/54klLSmO5MTYUzdTa3/a12WTc1hVuoemk62MyEzi0/OLGJeXPrAd08WF+N04E50m\nKTIIFOQkc9eNF1F2oJYX39zF1vJqAOobXSx7/G2q6xqpPd7EmcZVcTFuUj1x5GYkdgrsTiHuiSM1\nMTak2ze7XC5umJPP3Mk5vLB6J2u2VPDgM5u4qCiDG+cXkZ2eEOp/ukQgBbzIACjKTeHrn76YHftq\n+PdnNtHc0sq+quOkJsVSmJvSFtyBUXb760CIx52XufE0TxxLrivmymnDeea1MjaVHea9XUeYNyWX\nG2bnh/0WEtI/FPAiA2jMiFRSEmNxRcEPvxh5D4PJy07mn2++mI07fPzur2W8un4/a7dUcMOcfOZd\nnKt75Q9yCniR8yDKFf6HwfR0iwKXy8VUk8mkwgxe27CfFSXl/PbVD3h94wFumlekq3YHMQW8iACB\nh8Z/dMZILp2YzYtv7Wb1Owf56fObGTcqjU9fOZoRmUnhLlF6SX9/iUiQ5IRYbltg+PaSS5hYMJT3\n91Tzrd/8ncdeep/aY43hLk96QSN4EelWbkYiX7txMlt2HeGZ18t4491DlL5fxXWzRrFg+ojTnn8r\nkUcjeBE5owkFQ/n2HdO57WpDjDuK51fv4t5flvL39yuJlOtopHu60Kkf6eKNYOqPDk7pi/qGZlau\nLefV9ftobvFTmJvMp68cTWFOSq8+xyn90V8i6kInY0w08ChQ2PYZ37DWvtWlzfeAKwj8lfC/1tof\nhbIvEYkcCfHR3DiviCsuyuG5v+1kg/XxvSc2MLM4i09cXsjQlPhwlyidhDpFcxtw3Fo7B1gCPNR5\npTFmAjDPWjsbmA38H2NMdp8qFZGIkZmWwFc+PpF7bpnCqGwP67ZVcu+v1vHCGztpaOr53jhyfoUa\n8E8Bd7W99gFDu6yvBeKNMXFAPNAK1Ie4LxGJUGNGpHL/7dNYcu04EuOjWVmyh2/+Yh1vvnuQ1tZB\nP+s66PV5Dt4Y829Ai7X2/i7Lvwn8P8ANfMda+59n+pzm5hZ/tI7KiwxaDY3N/O/fynj+b2U0NrVQ\nkJPCkoXjmVTkDWq35LuvAPDofQvCUaYThT4Hb4xZCiztsvgBa+3LxpivAFOA67tsUwB8HCgAYoAS\nY8yz1tqqnvZTXT34B/g6cBRM/dHhQumLq6bkMnV0Bs+v3knJlgr+5b9LuHh0BjfOKyKr7UZmLS1+\n3G5X2PvDKQ8h93o9Pa47a8Bba5cDy7suN8YsIRDsi6y1Xe9rOh0otdbWt7XdDEwAXj/3skVkMErz\nxLH0umKunDqcZ1/7gHc+OMzmnUeYP2U4N8zJC3d5F5RQz6IpAO4ELrfWNnTTpAz4qjEmisAUzURg\nV8hVisigkz8smX++ZQobbOBGZn9Zv4+SLYfwE7iXvQy8UHt5KYEDq6uMMe3LFhA48LraWrvWGPMK\n0H7q5HJrbXlfChWRwcflcjFtbCaTizJ4dcM+VpaUc6KxhfqGZv7hx28QF+smPtZNXEzgZ3xsNHGd\n3p/62eV1fExbu1g38TGBn7HRUbopWhe60KkfXSjzrOdK/dFBfRFwtL6Je36+luYWP9npQ2hoaqHx\nZAuNTS00Nbf26bNdEBT67b8s2n9RdP5lEB/r5qXSvbhdLr73hZlhv/99RF3oJCISiuSEWBLjY3C7\nXXxnyYygda2tfhpPtgSFfkNT86llDU1ty9rWBV43d7Ms0LbmeBONTS1nren//uRNstMTKMxJpiA3\nhcKcZHK9ibijBv+dXBTwIhIRoqJcDImL7tcnWrX6/Zw82UpDU/OpXwLtv0B+tWIbLS2t5Ocks+vg\nUdZsqWDNlgog8NjE/GEeCnJSTgV/SmJsv9V1vijgRcSxolyuU9M2Xe+WExfjhhg33/j0xbS2+jl0\n5Dg7Dx5l54Fadh08it1bw/a9NafaZ6TEU5ibQkFOMoU5KYzMSuqXJ17d/UgJbrfr1EPZ+5MCXkQu\neFFRLnK9SeR6k7hscg4QuLHa7oqOwN95oJbSbZWUbqsEINodRV62JxD4bVM76cmRdS8eBbyInFcP\nfvnSQXHQOSE+mvF56YzPSwfA7/dTVX2CsvbAPxj4WXagFt7eBwSuAWgf4RfmJjMqy0NsTPiu0FfA\ni4icA5fLRVZ6AlnpCcyeOAyAxqYWyiuOnpra2XnwKBusjw3WB4A7ysWIzKRTgV+Qm4I3Jf68nc6p\ngBeRC1J/3KIgLtaNGZmGGZkGBEb5R442sPNAxwh/T0Ud5RV1vLYxsI0nIaYj8HNS2h6aMjCBr4AX\nEeknLpeLjJQhZKQMYUZxFgAnm1vYU3mMXW0j/J0Ha9lUdphNZYdPbTdQZ+go4EVEBlBMtJui3BSK\ncjvO46mua2TXwUDgv7ZhPwN1yv3gP5NfRGSQSfPEMdVkcuO8IpITYgfsAeYKeBERh1LAi4g4lAJe\nRMShFPAiIg6ls2hERMJoIK/s1QheRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdS\nwIuIOJQCXkTEoRTwIiIOpYAXEXGokO5FY4zJBB4H4oFY4C5rbWmXNrcAXwVagV9aax/tY60iItIL\noY7gbwWetNbOA+4FlnVeaYxJBP4VuAq4AviaMSa9D3WKiEgvhTSCt9Y+1OntCGB/lyYzgLettbUA\nxpg1wGxgRSj7ExGR3gv5dsHGmGwCge0B5ndZnQ34Or2vAoad6fPS0hKIHqDnEp5PXq8n3CVEFPVH\nB/VFMPVHsIHoj7MGvDFmKbC0y+IHrLUvA9ONMdcAjwELzvAxrrPtp7q6/mxNIt5A3dN5sFJ/dFBf\nBFN/BOtLf5zpF8NZA95auxxY3nmZMeZyY0yatbbaWrvKGPNEl80OEhjFt8sF1p17ySIi0lehHmRd\nDNwOYIyZCOzrsr6UwOg+1RiTRGD+/c2QqxQRkV4LdQ5+GfC4MWYxEAd8CcAYcw+w2lq7tu31y4Af\n+Hb7AVcRETk/XH6/P9w1AODz1UVGIX2gecVg6o8O6otg6o9gfZyD7/EYp65kFRFxKAW8iIhDKeBF\nRBxKAS8i4lAKeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdS\nwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIi\nDqWAFxFxKAW8iIhDRYeykTEmE3gciAdigbustaVd2twEfB1oBV6z1v5LH2sVEZFeCHUEfyvwpLV2\nHnAvsKzzSmNMAvBD4EpgFnCVMaa4L4WKiEjvhDSCt9Y+1OntCGB/l/X1xpiJ1to6AGPMEWBoyFWK\niEivufx+f0gbGmOygRWAB5hvrT3YQ7uJwLPAZGvtyZ4+r7m5xR8d7Q6pFhGRC5irpxVnHcEbY5YC\nS7ssfsBa+zIw3RhzDfAYsKCbbUcDTwM3nyncAaqr689WSsTzej34fHXhLiNiqD86qC+CqT+C9aU/\nvF5Pj+vOGvDW2uXA8s7LjDGXG2PSrLXV1tpVxpgnum5njBkO/AG4zVq7qfdli4hIX4R6kHUxcDuc\nmoLZ102bR4EvWWs3hrgPERHpg5AOshI4a+ZxY8xiIA74EoAx5h5gNXAEmAt8xxjTvs1D1to/9q1c\nERE5V6GeRXMYuLab5T/o9DYh1KJERKTvdCWriIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lAJeRMSh\nFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyI\niEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBwqOpSNjDGZ\nwONAPBAL3GWtLe2h7W+BRmvt50ItUkREei/UEfytwJPW2nnAvcCy7hoZYz4CFIa4DxER6YOQRvDW\n2oc6vR0B7O/axhgTB9wHfBdYHFJ1IiISspACHsAYkw2sADzA/G6afBP4b+BoqPsQEZHQnTXgjTFL\ngaVdFj9grX0ZmG6MuQZ4DFjQaZvRwDRr7beMMVecSyFpaQlER7vPte6I5fV6wl1CRFF/dFBfBFN/\nBBuI/nD5/f5eb2SMuRzYbK2tbnt/2Fqb0Wn9V4E7gHogGfACD1prf9TTZ/p8db0vJMJ4vR58vrpw\nlxEx1B8d1BfB1B/B+tIfXq/H1dO6UKdoFgMXAw8bYyYC+zqvtNY+DDwM0DaC/9yZwl1ERPpfqGfR\nLAM+Yox5A1gOfAnAGHOPMWZWfxUnIiKhC2mKZiBoisZ51B8d1BfB1B/BBmqKRleyiog4lAJeRMSh\nFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyI\niEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAK\neBERh1LAi4g4lAJeRMShFPAiIg4VHcpGxphM4HEgHogF7rLWlnZpMxl4tO3ti9baZX0pVEREeifU\nEfytwJPW2nnAvUB34f1L4AvAJUCxMSYhxH2JiEgIQhrBW2sf6vR2BLC/83pjTBaQZK3d2LboM6GV\nJyIioXL5/f6QNjTGZAMrAA8w31p7sNO6GcC/AzuB0cBz1tqHz/R5zc0t/uhod0i1iIhcwFw9rTjr\nCN4YsxRY2mXxA9bal4HpxphrgMeABV12mA8sAk4Aa40xf7HWbu1pP9XV9WcrJeJ5vR58vrpwlxEx\n1B8d1BfB1B/B+tIfXq+nx3VnDXhr7XJgeedlxpjLjTFp1tpqa+0qY8wTXTarBLZaa4+0tX8LGA/0\nGPAiItK/Qj3Iuhi4HcAYMxHY13mltXY34DHGpBtjooCLANuXQkVEpHdCOshK4KyZx40xi4E44EsA\nxph7gNXW2rXA14CXAD/wZ2vtu/1Qr4iInKOQD7L2N5+vLjIK6QPNKwZTf3RQXwRTfwTr4xx8jwdZ\ndSWriIhDKeBFRBwq1Dn4fnf3IyXdLv/ojJFcOXU4AL9asZUd+2pPa1OYm8ydCycAsHrTAVaW7On2\ns77/xZlEu6M4dOQ4Dz3b/SGBz10zlvF56QAse/xtjh4/eVqb2ROzWTS3AIBnX/+A9dt9ALjdLlpa\nAjNN3tR4/unmKQC8s8PH069+0O3+vnnrFNKT4znecJJv/frtbtt84vICZo7PBuDh597lgO/4aW0m\nFqTz2Y+OBeCldXt4feOB09rExbr57tIZAOw8UMvPX+z+pKavLJ5AXnYyAPf8Yu2p/6bOPjJ9BAum\njwDg16ve5/3y6tPajBmVxuevHQfAW5sP8eJbu7vd37KllxAfG01VzQkefPqdbtvcdrVhUuFQAP7t\nqQ1UH208rc2M4iw+eUUhAL//205Kt1We1iYtOY57b50KwOadR3jy5e6P/d9988Vkpg6hoamZ+5f/\nvds2C+fkM2fSMAB+9sJ77Kk4/U/scXlp3HFNoA9eeXsff3l732lt3G4XP/jiLADKK47yXy9s6XZ/\ndy4cT2FuCgD3LS+lsanltDbzp+TysZmjAHjiz9t5b9eHp7XJ9Sby1U9NBmDd1gqeX72r2/19647p\nJMbH8OHRBr7/1MZu29x81WguHuMF4EdPb8RX03Bam2ljvdw0fzQAf3hzF2veqwj6/wpAcmIM998+\nHYCt5R/y2Krt3e7vrpsmM2xoIs0trXzzF+u6bXPdpaO4/KJcAH7+4hZ2Hjh6WpsxI1L4/PXjAXht\nw37+XLq328968MuXArC/6hg/+f3mbtssvW4cZmQaAA/8+u/UNzSf1uayycO4fnY+AP/zlx1s+uBw\n0Hq324U3dQhfv+kiANZvr+LZ18u63d99t08jJTGW2uNNfPfx9Tz2wNXdtgON4EVEHEsHWfuRDhwF\nU390UF8EU38E00FWERHpFQW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lAJeRMShFPAiIg4VMRc6\niYhI/9IIXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4ERGHiphH9g12xpgfAXMJ9On3\nrbUvhLmksDLGDAG2AMustY+FuZywMsbcAvwT0Az8q7X2T2EuKWyMMUnAE0AaEAd821r7cnirOv+M\nMROAF4EfW2t/ZowZATwJuIFDwG3W2tOfS9lLGsH3A2PMPGCCtXYW8FHg4TCXFAnuA05/IOgFxhgz\nFHgAmANcBywMb0Vh9znAWmvnAZ8EfhLecs4/Y0wi8J/Aa50Wfwf4L2vtXKAMuKM/9qWA7x9vAJ9q\ne10DJBpj3GGsJ6yMMWOBYuCCHal2chXwqrW2zlp7yFr7hXAXFGaHgaFtr9Pa3l9oGoFrgIOdll0B\n/LHt9QoC35s+U8D3A2tti7X2eNvbJcAqa+3pj7y/cPwHcFe4i4gQeUCCMeaPxpg3jTFXhrugcLLW\nPgOMNMaUERgYfSPMJZ131tpma+2JLosTO03JVAHD+mNfCvh+ZIxZSCDg/yHctYSLMeazwFpr7e5w\n1xIhXARGrIsJTE/8xhjT40OSnc4Ycyuw11pbBMwHfhbmkiJRv30/FPD9xBhzNfAvwMestbXhrieM\nrgUWGmPWAUuB+40x/fLn5iBVCZS0jdp2AnWAN8w1hdNs4GUAa+27QM6FPJ3ZybG2ExMAcgmevgmZ\nzqLpB8aYFOBB4Cpr7QV9YNFae1P7a2PMt4Bya+2r4aso7F4BHjPG/JDAnHMSF+a8c7syYAbwvDFm\nFHDsAp/ObPcq8Angqbaff+6mHxg7AAAAhklEQVSPD1XA94+bgAzgd8aY9mWftdbuDV9JEgmstQeM\nMb8H1rUt+kdrbWs4awqzXwC/NsasJpA/d4a5nvPOGDOVwHGqPOCkMeaTwC0EBgJfBPYAj/fHvnQ/\neBERh9IcvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIO9f8B84EwBdn/DFIAAAAA\nSUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f26b8b2c908>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "Binning Statistics: Energy = -2.9215 +- 0.0288, Auto correlation Time = 0.4686\n",
            "Step 10, dE/|E| = 0.1998, elapse = 2.4429\n",
            "Binning Statistics: Energy = -3.0027 +- 0.0277, Auto correlation Time = 0.4469\n",
            "Step 11, dE/|E| = 0.1776, elapse = 2.4277\n",
            "Binning Statistics: Energy = -2.9655 +- 0.0311, Auto correlation Time = 0.5804\n",
            "Step 12, dE/|E| = 0.1878, elapse = 2.4272\n",
            "Binning Statistics: Energy = -3.0553 +- 0.0278, Auto correlation Time = 0.5056\n",
            "Step 13, dE/|E| = 0.1632, elapse = 2.4185\n",
            "Binning Statistics: Energy = -3.0617 +- 0.0290, Auto correlation Time = 0.5580\n",
            "Step 14, dE/|E| = 0.1614, elapse = 2.4224\n",
            "Binning Statistics: Energy = -3.1071 +- 0.0278, Auto correlation Time = 0.5621\n",
            "Step 15, dE/|E| = 0.1490, elapse = 2.4376\n",
            "Binning Statistics: Energy = -3.0956 +- 0.0242, Auto correlation Time = 0.3862\n",
            "Step 16, dE/|E| = 0.1521, elapse = 2.4204\n",
            "Binning Statistics: Energy = -3.1836 +- 0.0257, Auto correlation Time = 0.4730\n",
            "Step 17, dE/|E| = 0.1280, elapse = 2.4350\n",
            "Binning Statistics: Energy = -3.2158 +- 0.0231, Auto correlation Time = 0.3985\n",
            "Step 18, dE/|E| = 0.1192, elapse = 2.4409\n",
            "Binning Statistics: Energy = -3.1611 +- 0.0270, Auto correlation Time = 0.5128\n",
            "Step 19, dE/|E| = 0.1342, elapse = 2.4217\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XeYVNX9x/H3lO19YWHpnaPs0kGK\nBbFgj4oae0wiiS2JDRL1l0QNxpCQGBOjMQSMir0LFrCj9N7hSO9ll7Kwu2yf3x8z4LjMwjK7O7vc\n/byex4eZe86d+51h/HDm3Oby+XyIiIjzuOu7ABERqRsKeBERh1LAi4g4lAJeRMShFPAiIg7lre8C\nDsvJOdhgD+dJS4tn377C+i6jSqqvZhp6fdDwa1R9NVOT+jIyklxVtWkEXw1er6e+Szgm1VczDb0+\naPg1qr6aqav6FPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoU76gB/1\nzExGPTOzvssQEWlwTvqAFxGR0BTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUCd9wJeVV1BYVEpB\nUWl9lyIi0qCc9AFfWlbBoZJyHv3fPNZvPxDx7es4fBFpqMK6ZZ8xxgtMADoFXmOktXZ6pT49A30A\n3rfWjq5JoVWJjfbg8/nYk1fEn15awDVDO3N+v9a4XFXexUpEpFEIdwR/M1BgrT0DuBV4IkSfccDP\ngdOAbsaY+DC3dUwul4v42Cjuu64XCXFRvPb5Gv71zjJN2YhIoxfuTbdfAl4NPM4BmgQ3GmOaA4nW\n2oWBRdeHuZ1qy2qfzqM/6c+4yStZtCaXR56bx+1XZNGpZUpdb1pEpEEKawRvrS211hYFnt4DvFKp\nS3tgrzHmeWPMDGPMPTWosdpSEmO4/9peXH5GB/YeKGLMSwuZOnczPp8vEpsXEWlQjjuCN8aMAEZU\nWvywtXaqMeYuoA9wWaV2F9ABuAI4BMwyxnxqrV1R1XbS0uLDurO4x+Ofa8/ISDqybMSVPeif3YK/\nvryA179Yy8Zd+dx9XW+S4qOPWv/Wxz4BYMJvhx1zO8Gvf7ztR1J1669v9fX5VFdDrw8afo2qr2bq\nor7jBry1djwwvvJyY8yt+IP9Cmtt5QnvXcAKa+2eQN/pQBZQZcDv21d4AmV/p7zcPzrPyTn4veUt\nU2N5+JZ+jJu8kjkrdvLLsV9w2+XZdG6VUq31g2VkJFXZXp3161J5uQ+Px1Vv26+OY31+DUFDrw8a\nfo2qr2ZqUt+x/mEIa4rGGNMRuB0YHjRVc4S1dgOQZIxJN8a4gV6ADWdbNXF4yuaKMzuw92Axf355\nIVPmbKaihlM2xaXlLFu/h9c+X8P+/GL2Hiziq0XbNBUkIg1KuDtZR+DfsfqRMebwsmHAfcA0a+0s\n4F7gY8AHTLHWLqlhrSGNvXPwMdvdbhc/OL0DXVun8p9JK3jjy7XYzfu49dJuJMZFVWsbFRU+Nu86\nyIoNe1m+YS9rtu6nrPy7MHcBL061LF6by48vOoXUxJiavCURkVoRVsBbax8CHgrRNCaozxxgQJh1\n1bpT2qXxyE9PY/zkFSxZt4dH/jeX23+QXWX/vPxiVmzcy4oNe1m1eT/7DxYfaWvbLJGsDulkd0hn\nwoerqPBBq6bxLF23h99PmMuPLjD0O6VZJN6WiEiVwh3Bn5RSEqK599pefDhrE+99s54xLy8kNtpD\nbLSH0rJyvt2Sd2SUvjUn/8h6aUkxDM7OJKtDOt3ap5OS8N3OWpfLhccF917biy8XbuONL9fyzHvL\nGZydyQ3ndSU+tlF9xCLSgDS69HG7XFw2uD1dW6fw7KQV5OWXUFxazi+e/IbSsgoAvB43We3TyOrQ\nhOwO6fTqlklubv5xX/fcvq3p1j6N/05eyczlO/1TQZd045R2aXXyXvbnF1NYVEq5z0fO/kNkpMbV\nyXZE5OTU6AL+MNM2jUd/chqjnplJaXkFLZrEk9UhnawO6XRtnUp01HeHbJ7IZQ9aNEngoZv78sHM\njXwwcxNjX13E+f3bcNWQjkSFOAz08HVsjrcvIdi67Xl8Pn8r81bvprzCvy/gr68t4sGb+mr+X0SO\naLQBD5CcEE1yQjQ+n48/3Fp7uwu8HjdXnNmRHp2a8t/JK/hk3hZWbNjLiEu70S4zvGNdy8ormLd6\nN5/N38qGHf6LqrVoEs+BghJ8QM7+Iv72+mJ+c0Ofau88FhFnO+mvJlkb6urCZB1bJvPIT05jaJ9W\nbMst4LEX5/PhrI1UVFT/cMq8/GLen76BUc/M5L+TV7JxxwF6dW7KyOt68diIAcRGe0mI9XJu39Zs\nyyngyTeXUFRSVifvR0ROLo16BB8JMdEebh5m6N25KRM+WsXb09azZN0eRlzajWbHmDNfv/0Any3Y\nwrxV/mmYuBgvw/q34Zy+rY9az+Vycf15XThUXMbM5Tt56u1l3HNNj5BTQiLSeCjga6i6c+fZHZsw\n+tYBvDjVMn/1bh6eMJfrzu2Mz+c78guirLyC+at389mCrUeubd+iSTzn9WvDoKzmxEZX/dfldrn4\nycWncKi4jEVrcnn2/RXceWU2Hrd+pIk0Vgr4CEqMi+KOy7OY3aUpL33yLS9MsUR53cTHeJk0fQNf\nLtpGXkEJLqBX56ac26813dqlVXsKyeN2c/vlWTz55lIWrcnlfx+t5qeXnIpb18YXaZQU8BHmcrkY\nlJWJaZPKhA9XsWrTPvLKSnhv+obvpmH6tKJZWniXz4/yevjF8O787fXFzFy+k7gYLzec10U3QBFp\nhBTw9SQ9OZb7r+vFr578hpKycq4/twuDsjOPOQ0Tytg7Bx91oaK4GC/3XNOTP7+ykM8XbCUh1ssV\nZ3as7bcgIg1cow/4Ezn+vLa5XS7iYrzExXgZ2qd1rb52YlwU91/biz+9tIBJMzYSH+Nl2Glta3Ub\nEN5x/CISGdoD52CpiTGMvK43qYnRvPbFWr5Zsr2+SxKRCFLAO1xGahz3X9ebxLgonp+ymvmrd9d3\nSSISIQr4RqBV0wTu/WFPoqM8/GfSCpav33OkbdQzM49Ms4iIsyjgG4kOLZK5+6oeuFwu/vXuMtZu\nzauV1/X5fDW+gYqI1A0FfCNySrs07rwim7IyH39/cwmbd534LcLyCkpYtCaHt6etY+yri9h3sJh9\nB4uZu2pXHVQsIjXR6I+iqW+RPvqkV5emjLj0VP47eSVPvL4Yt8uFxxP63/nSsnI27cxn/fY81u84\nwPrtB8jN+/4dGt1uF74KHxM+XEXTlDg6tkyOxNsQkWpQwDdCA7MyOVRcxsRPvsXtguSEGHw+H7v2\nHfKH+XZ/mG/ZnX/kcsTgP/SyR6cmdGyZTMeWyXRokcwjz82jpLSc/KJSnnp7Kb+7pR/pybH1+O5E\n5DAFfCM1tE9rCovLeHvaevIKivnVP76hoOi7q1B6PS7aNk+iUyDMO7ZMJiM1LuQZsdFRHq4b3J5X\nP1/DP99ayoM39SUmWhc6E6lvCvhG7OKB7fhw1iaKSsqJj/WS3fG70XnbZklEeau/i+a8fq3ZsaeA\nrxZvZ9zkFdw1vLuugSNSzxTwjZjL5SIhNor4GC9/vr1m+wJcLhc3nN+VXfsOsWhNLu9+vZ6rhnSq\npUpFJBwKeKnRhciCdxJ7PW7uuCKbP744nw9nbSIzPZ7Tu7eojRJFJAw6TFJqVWJcFL+6ugfxMV5e\nmLKaNVv313dJIo2WAl5qXYsmCdxxZTYVFfCvd5axa29hfZck0igp4KVOZLVP58ZhXTlYWMroCbM5\nVKz7xIpEmgJe6szQ3q04t29rNu08yH8mrTihm42LSM1pJ2sjV9dn0l53bmf25Zew0O7mjS/Xct25\nXep0eyLynbAC3hjjBSYAnQKvMdJaO71Snz8CZ+P/lfCutfYvNStVTkYet5tf39yPe//+FZ/M20LL\npgmc1bNlfZcl0iiEO0VzM1BgrT0DuBV4IrjRGJMNDLXWng6cDvzEGJNZo0rlpJUQF8XdV/cgMS6K\niVMtqzftq++SRBqFcAP+JeC+wOMcoEml9jwg1hgTA8QCFYAOpWjEmqXFc9eV2QA8/e4ydu3T10Gk\nrrl8NbyWtzHmcaDcWvu7SssfBO4GPMAfrLVPHet1ysrKfV6vrl/idJ/O2cQ/31hM62aJjP3VWdz9\nty8BmPDbYfVcmchJq8ozFY87B2+MGQGMqLT4YWvtVGPMXUAf4LJK63QErgQ6AlHATGPM69baKu8X\nt68Bj+gyMpLIyTnxa6dHyslUX6+O6Vx4WlumzN3MY+NnUVZWgcvlCrv+2rjpd0P//KDh16j6aqYm\n9WVkJFXZdtyAt9aOB8ZXXm6MuRV/sF9hrS2t1NwfmGOtLQz0XQpkA19Uv2xxqqvP7sSOPQUsWbeH\nmCgPiXFR9V2SiCOFNQcfGKHfDgy31haF6LIW6GeMcRtjooDuwPrwyxQncbtd/PwHWbTOSKC4tJyi\nEp0EJVIXwj0OfgT+HasfGWMOLxuGf8frNGvtLGPMJ8DhQyfHW2s31qRQcZa4GC+/uroHv3l2FgVF\nZdz71HRSE2NISYwmNTGG1MCf3z2PITkhCo9b5+aJVFdYAW+tfQh4KETTmKA+DwMPh1mXNAJNU+JI\njo+moKiMmCgPO/YUsOkY94l1AUkJ0UfCPzUxmsKiUtxuF4VFpcTHaqpHJJjOZJV65fW4SUmIZszt\ng/D5fBwqLmd/fjF5+cXszy9hf0Ex+w+WkFcQeJ5fzM69hWzelf+91xn5zEzO7tWK8/u3IS0ppp7e\njUjDooCXBsPlchEf6yU+1kvLpglV9vP5fBSV+P8h+PPLCykr9xEV5WbK3M18On8Lg7IyuXBA22O+\nhkhjoICXk47L5SIuxktcjJcor4coLzz+84HMWrGTKXM2M33ZDqYv20HvLk25aGA7OrdKqfK1Rj0z\nE4/HxZjbBkXwHYhEhgJeHCHK6+asni05o0cLFq/J5ePZm1i0JpdFa3Lp0jqFiwa2o0enJrpPrDQq\nCnipV7V9NUu3y0Wfrhn07tKUb7fs5+M5m1m6bg9r3lpKq6YJXDigLQO6Ncfr0dE44nwKeHEkl8uF\naZuGaZvG1t35fDxnM3NX7WLCh6t45+v1XNC/DWf10lUtxdkU8HJSq84vgNbNEvnZZd0YflZHPpm3\nha+XbOe1L9YyeeZGKip8JOhMWnEo/U6VRqNJSizXn9eFsXcO5oozO+ByuThUUk5uXhEvTrXsbsDX\nQxIJhwJeGp3EuCh+cHoHxt45mIRYL263i68WbePBcbN55r3lbNhxoL5LFKkVmqKRRismykNstP+4\n+6uGdOLj2ZuZv3o381fv5tR2aVw0oC1ZHdJx6cgbOUkp4KXRc7lcnHZqc/qf0oxVm/bx8ZzNrNiw\nl1Wb9tGmWSIXDWhL/1ObhbwOTm1crlikrijgRQJcLhfd2qfTrX06m3YeZMpc/5E34yav5O1p6xl2\nWhvO6tGSmGjdmEZODpqDFwmhXWYSt/0gizG3DeLcPq05WFjCq5+tYeQzM3jvm/UcKCyp7xJFjksj\neJFjyEiN48ZhXfnBGe35YuE2Pl+wlUkzNvLxnM2c0aMF5RUVuoSxNFgKeGnUxt45uFq3S0uKj+by\nMzpw4YC2TF+6g6lzN/Plwm0ARHvd5OUXk5Koq1hKw6Khh8gJiInycG7f1vzptoHc9oMsPG4XJWUV\n/HHiAnbu1XH00rAo4EXC4HG7GdCtOSkJ0cTFeMnNK+LxiQtYty2vvksTOUIBL1IDLpeL+BgvP77o\nFAqLyhj76iIWr82t77JEAAW8SK04q2dLfnlVdwCeenspXy/ZHrFtj3pm5pHj8UWCKeBFaknPzk0Z\ndUNvEmKjeP7j1UyavgGfz1ffZUkjpqNoRGqg8hmsnVqm8NDNfXni9cW8N30Dew8Wc/MFXXUopdQL\nfetEallmejz/d3Nf2jZP5Osl23n6neUUl5bXd1nSCCngRepASmIMv7mhD1nt01i8Npe/vraI/EOl\ntbqNvIISPpy1kbyCYkrLKmr1tcUZFPAidSQuxsvd1/RkYFZz1m07wOMTF5C7/1CNXtPn82E37+PZ\n95cz8ukZvD1tPWXlPg4UlrBq495aqlycQgEvUoe8HjcjLu3GRQPasnNvIX+cuIDNu4591mwoBYdK\n+XzBVn43YS5/fmURc1ftpnl6PDee35XEwB2pnnxrKSsU8hJEO1lF6pjb5eKaoZ1JTYrhtc/WMObl\nhfxyeHdObZ8OHPuSw5t2HuTLRduYs2oXxSXleNwuTju1GUN7t6Jrm1RcLhdT5mzGBRwqKeefby3l\nl1d1J7tDk0i+RWmgwgp4Y0wz4AUgFogG7rPWzqnU50bgHqACGGetnVDDWkVOauf3a0NKQjTjP1jJ\nE28s4dZLT2Vgt8yj+pWUljNv9W6+XLSN9dv9d5dqlhbHmYNacEaPlqQkRB+1TnSUh59d1o2n3lnG\nP99axi+v6k73jgr5xi7cEfxNwERr7SvGmCHAaGDY4UZjTALwe+A0oASYZ4x511qr34/SqJ12anOS\n46N56p1ljJu0kv0Hv7vs8K69hXy5aBszlu2goKgMF9CjUxOG9m7F0AHt2bsn/5ivnd2xCb+6ugf/\nfGspT729lLuu7E7Pzk3r+B1JQxZWwFtrnwh62gbYWqnLAGCetTYPwBgzAzgdmBzO9kSc5JR2aTx4\nYx+eeGMxb3y5lpgoDxUVPh4cNxuA5PgoLhnUjiE9W9I0NQ4Aj7t6tw3Map/OPVf34B9vLeVf7yzj\nriu706uLQr6xcoV7pp0xJhN/YCcB51hrtwe13QD0t9beG3g+GthirR1X1euVlZX7vF7dKUcaj937\nCnnkv7PYsss/Ms/q2ISLB7dnUPeWRHlrdvzDsrW5PDphNuXlFfzmR/0ZmN2iNkqWhqnKf/2PO4I3\nxowARlRa/LC1dirQ3xhzMfA8QVM0J1LAYfv2NdxLrVbneuH1SfXVTH3V5wJ+fX1vfvPvWUR53dz/\nw54A7N9XcFTfE60xMyWGe67uwZNvLmXMC/O4/fJs+pqM2iq9xvVFmpPry8hIqrLtuMMEa+14a+3A\n4P+AImNMWqD9I6BPpdW2A8F7j1oFlolIkITYKOJivHg9tX/Esmmbxr0/7InX6+bZ95czf/XuWt+G\nNGzhfquGA7cAGGO6A1sqtc/BP7pPNcYk4p9//ybsKkUkLF3bpHLfD3sS5XXz7PsrmLtqV32XJBEU\nbsCPBs43xnwNjAfuADDGPGCMGWStPQQ8AEwFPgMePbzDVUQiq0vrVO67thfRUW7GTVrJnJUK+cYi\n3KNocoFLQiwfE/T4LeCt8EsTkdrSuVUK91/XiydeX8y4ySuo8PkYlOWfRT3WiVZyctOZrCL1LFLB\n2qllCiOv681fX1vM+A9W4vP5GFwLR9eMemYmHo+LMbcNqoUqpTbpWjQijUiHFsmMvK4XcdFeJnyw\nihnLdtR3SVKHFPAijUyHFsmMur438bFenvtwFUUlZfVdktQRTdGINELtMpMYdb1/uib/UCmFxWU8\nNG42CbFeEuKi/H/GRlV6/P1l8bFe3amqgVPAizRSbZv7Q/6xF+ZT4fNRWFRKzv5DlFdU/+z2uBgP\nxaUVxEV7KCwqJT42qg4rlhOlgBdpxNo0SyQ5cHXKsXcOxufzUVxaTsGhMgqKSik4VEpBUeBxUVng\nefDjMrbuzqegqIxR/57FsP5tOL9fawV9A6GAF5EjXC4XsdFeYqO9NEmJrdY6I5+eQXFZOS5cvD99\nA5/M26KgbyA0gSYiNeJyuUiIjeIvdwzimrM74XH7g37Uv2fx/vQNFBbV7r1oKxv1zMwjx/LL92kE\nLyK1Ijbay0UD2zG0Tyu+XLiNj+dsrtaIXida1R0FvEgjV9vBGm7QS+1TwItIjYy9c3DIy93WZtBX\nVPjYn19Mbl4RuXmH/H/u9z/en19MXIyiLBR9KiJSp44X9BU+Hy5g38Hi78I7r4jc/f7He/KK2HOg\n6JiHb+YfKuWbJds5s2fLyL2xk4ACXkQioqqgP+z+p2eEXC85IZp2mUk0TYmlaUpc4M9YmqbG0SQ5\nhl//exYHCkv438erKSuvYGif1pF6Sw2eAl5EIqpy0L/z9XoA+nTNOBLeTVLiyEiNJT05lpioY9/K\n0+txkxwfjc/nY+In31JaVsGw09pG4q00eAp4EakXh4P+i4XbALjjiuywX8vrcXPftT0Z++oiXvti\nLaXlFVwyqH0tVXryUsCLyEkt+CigB27sw9hXF/H2tPWUllVw+RkdcLmOe0tox9KJTiLiGM3S4vnN\njX3ISI1l0oyNvPXVOny+6l9bx2kU8CLiKE1T4njgxr40T4/n4zmbefWzNY025DVFIyL1qi7OYE1L\niuGBG/yXQ/5swVa80V6uPqsD7gY4XVOXd8TSCF5EHCklMYZf39Cbts0SmTJrI//7aBUVJ3ApZCdQ\nwIuIYyXFRzPqht50aZPKjGU7GTd5BWXlFfVdVsQo4EXE0RJio3js9sF0bp3C3FW7efb9hhPyFT4f\npWUVdbaPQAEvIo4XHxvFfT/sySltU1n4bQ7/emcZpWXlQP1dbnjjzgM8PnEBBwpLOFRcXifbUMCL\nSKMQG+3lnmt6kt0hnaXr9vDPt5ZSXFo3wXos+YdKmTjVMvr5+azffoBor5vY6GOfrRsuBbyINBrR\nUR5+eVUPenVuyoqN+/j7G0sidghlhc/HN0u289C42Xy5aBuZTeIZdV0vkuKjcbvr5ugeHSYpIo1K\nlNfNnVdmM27SCubbHLweF4lx0XW6zU07D/LSp5Z12w4QE+XhmqGdOL9fG7yeuh1jhxXwxphmwAtA\nLBAN3GetnVOpz7XA/UAF8Lm19v9qWKuISK3wetzcdnkW3g9XMXvFLvbnFzPm5YUMzGpOP9OMxLja\nuSFJYVEp7369gS8WbcXng/6nNOPaczqTnly9+93WVLgj+JuAidbaV4wxQ4DRwLDDjcaYeODPQHcg\nH5htjHnZWruypgWLiNQGj9vNiEu6sWzdHopLy/l2y36+3bKflz/5lu4dmzAwqzk9Ozc97tUsQ/H5\nfMxcvpM3v1zLgcJSmqfHc9P5XcnqkH5U36pumFIbwgp4a+0TQU/bAFsrtRcaY7pbaw8CGGP2AE3C\nrlJEpA643S5io73ERnt58KY+zFm5i9krd7F4bS6L1+YSE+2hb9cMBmY159R2aXjcR0+pVL6n7Jbd\n+bz0iWXN1jyio9xcNaQjw/q3Jcob+V2eYc/BG2MygclAEnBO5fagcO8OtAdmH+v10tLi8XrrZk9y\nbcjISKrvEo5J9dVMQ68PGn6NJ2t9Ho9/B6fplIHplMGPLstm084DTFu4lWkLtzJz+U5mLt9JalIM\nZ/Zqxdl9WtOlTeqRq1QeXj8+MZZXpq7mgxkbqKjwMah7C0Zcnk2ztPga1VcTruPtQTbGjABGVFr8\nsLV2aqD9YuAea+2wEOt2Ad4BbrbWLj7WdnJyDjbYc4jr6udTbVF9NdPQ64OGX+PJXF/lEXiwCp+P\nddvymL1iF/NW7yb/UCkAzdLiGNitOQOzMvnrq4soKa3A7XFxoKCEZmlx3Hh+V7p3rP6kRU0+v4yM\npCoPwTluwIcSmHdfaq3dF3iea61tWqlPa2Aq/nBfeLzXVMCHT/XVTEOvDxp+jY2hvrLyCpZv2Muc\nlbtY9G0OJWX+s2FdLvD5/EfnXDqoHRcOaEvUCc5G1FXAhztFMxzoDTwZmILZEqLPBOCO6oS7iEhD\n5/W46dW5Kb06N6WopIxF3+Yye+Uulq3fQ5TXzR9HDKBpalx9l/k94Qb8aOAFY8xwIAa4A8AY8wAw\nDdgDnAn8wRhzeJ0nrLWTalauiEj9i432Mig7k0HZmYx8egYul6vBhTuEfxRNLnBJiOVjgp5Wb8+C\niMhJrCHfElCXKhARcShdqkBEpAbq4o5UtUUjeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiI\nQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4\nERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lDecFYyxjQDXgBigWjgPmvt\nnCr6vgoUW2t/HG6RIiJy4sIdwd8ETLTWDgUeAkaH6mSMOR/oFOY2RESkBsIawVtrnwh62gbYWrmP\nMSYG+C3wGDA8rOpERCRsYQU8gDEmE5gMJAHnhOjyIPBv4EC42xARkfC5fD7fMTsYY0YAIyotftha\nOzXQfjFwj7V2WNA6XYC/W2svNcacDfz4eHPwZWXlPq/Xc+LvQESkcXNV2XC8gA/FGDMEWGqt3Rd4\nnmutbRrUfg/wU6AQSAYygLHW2r9U9Zo5OQdPvJAIychIIifnYH2XUSXVVzMNvT5o+DWqvpqpSX0Z\nGUlVBny4UzTDgd7Ak8aY7sCW4EZr7ZPAkwBBI/gqw11ERGpfuEfRjAbON8Z8DYwH7gAwxjxgjBlU\nW8WJiEj4wj2KJhe4JMTyMSGWfQV8Fc52REQkfDqTVUTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoB\nLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4\nlAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAX\nEXEobzgrGWOaAS8AsUA0cJ+1dk6lPj2BCYGn71trR9ekUBEROTHhjuBvAiZaa4cCDwGhwnsc8HPg\nNKCbMSY+zG2JiEgYwhrBW2ufCHraBtga3G6MaQ4kWmsXBhZdH155IiISLpfP5wtrRWNMJjAZSALO\nsdZuD2obAPwVWAd0Ad601j55rNcrKyv3eb2esGoREWnEXFU1HHcEb4wZAYyotPhha+1UoL8x5mLg\neWBYpQ12AK4ADgGzjDGfWmtXVLWdffsKj1dKvcnISCIn52B9l1El1VczDb0+aPg1qr6aqUl9GRlJ\nVbYdN+CtteOB8cHLjDFDjDFp1tp91tqPjDEvVlptF7DCWrsn0H86kAVUGfAiIlK7wt3JOhy4BcAY\n0x3YEtxord0AJBlj0o0xbqB1Kk+bAAAIwElEQVQXYGtSqIiInJiwdrLiP2rmBWPMcCAGuAPAGPMA\nMM1aOwu4F/gY8AFTrLVLaqFeERGppnCPoskFLgmxfEzQ4znAgPBLExGRmtCZrCIiDqWAFxFxqHDn\n4GvdqGdmhlx+4YC2nNu3NQD/nbyCb7fkHdWnU6tkbr88G4Bpi7fxwcxNIV/rT7cNxOtxs2NPAU+8\nHnqXwI8vPoWs9ukAjH5hHgcKSvF4XJSXf3e+wOndM7nizI4AvP7FGuavzjnqdTJSY/n1DX0AWPRt\nDq98tibk9h68qQ/pybEUFJXyyHPzQva5akhHBmZlAvDkm0vYllPwvXaPx0W3dmn86MJTAPh49ia+\nWLjtqNeJifbw2Aj/rNm6bXk8+37og5ruGp5N+8xkAB74z6zvvffDzu/fhmH92wDw3EerWLVx31F9\n2mUm8Yvh3QGYvnQH70/fEHJ7o0ecRmy0l937DzH2lUUh+9x8gaFHpyYAPP7SAvYdKD6qz4Buzbn6\n7E4AvPXVOuas3HVUn7TkGB66qS8AS9ftYeJUe9TfL8CoG3rTLDWOopIyfjd+bsiaLj+jA2f0aAHA\nv95ZxqadRx/mdmr7NH568akAfDJvC5/O23JUH4/HxZjbBgGwcecBnn5necg+P7u0G51apQDw2/Fz\nKC4pP6rfOX1acdHAdgC8OGU1y9bvPapPq4wE7rmmJwCzV+zk7WnrQ76/R37an4TYKPYeKOJPLy0M\n2eeG87rQu2sGAH95ZSE5+4uO6tPvlAyuPacLAO99s54Zy3Ye1Sc5IYrf3dIfgBUb9/L8R6tDbu++\na3vSokkCZeUVPPif2SH7XDq4HUN6tQLg2feXs27bgaP+jru2SeFnl2UB8PmCrUyZsznka429czAA\nW3fn84+3lobsM+LSUzFt0wB4+Lm5FBaVHdXnrJ4tuOz0DgC8/Om3LF6T+712j8dFRmoc91/bC4D5\nq3fz+hdrQ27vt7f0IyUhmryCEh57YT7PP3xByH6gEbyIiGOFfSZrbcvJOdgwCgnBySdJRILqq7mG\nXqPqq5kanuhU5ZmsGsGLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh2ow\nJzqJiEjt0gheRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdqMLfsayiMMX8BzsT/\n2fzJWvtOUNtGYAtw+F5pN1prj743Xt3VdjbwJnD4XnvLrLW/DGo/D3g8UN9H1trRkaotsP1bgZuD\nFvWz1iYGtZcCM4Laz7XWHn3fudqvKxt4H/i7tfZfxpg2wETAA+wAbrbWFlda5+/AQMAH3G2tDX0/\nxbqr739AFFAK3GSt3RnU/2yO8T2IQH3PA32BPYEuY621H1ZaJ2KfXxU1vglkBJrTgdnW2p8H9f8x\nMBpYF1j0qbX2j3VU2/cyBZhHhL5/CvggxpihQLa1dpAxpgmwCHinUreLrLX5ka/uiGnW2quraPsn\ncAGwDZhmjHnbWrsyUoVZaycAEwCMMUOAH1bqkmetPTtS9QTqSACeAj4PWvwH4Glr7ZvGmMeBnwL/\nDlpnCNAl8D04FXgOGBTB+h4Dxllr3zDG3AXcB/y60qrH+h7UdX0AD1prP6hinYh9flXVaK29Jqj9\nOWB8iFVft9aOrKu6AtsOlSmfE6Hvn6Zovu9r4PAXYz+QYIzx1GM91WaM6QjstdZusdZWAB8B59Zj\nSb/HP0Kqb8XAxcD2oGVnA5MCjycD51Va51zgPQBr7SogzRiTHMH67gTeDjzOAZrU0barI1R9xxPJ\nzw+OUaMxxgCp1trQd06ve0dlChH8/mkEHyQwXVAQeHor/mmOylMIzxpj2gPT8Y9iIn2th27GmEn4\nf3Y+aq39NLA8E38YHLYb6BTh2gAwxvQHtgRPKwTEGmNeAdoBb1trn6jrWqy1ZUCZ///zIxKCfhLv\nBlpUWi0TWBD0PCew7EAk6rPWFgAEBhd34f/FUVlV34M6ry/gF8aY+/B/fr+w1uYGtUXs8ztOjQB3\n4x/dhzLEGDMF/1TYSGvtojqo7ahMAS6I1PdPI/gQjDGX4//L+EWlpt/j/7l8NpANXBXZylgDPApc\nDtwCTDDGRFfRt8ob8UbACOD5EMtHAj8HhgE3GmP6RbKoKlTnc4r4ZxkI94nAF9baytMjJ/I9qAsT\ngQestecAi4FHjtO/Xr6Lgc/kDGvtlyGaZwOPWGsvBH4LvFjHtVSVKXX6/dMIvhJjzAXA/wEXWmvz\ngtustS8G9fsI6A68FanaAjt0Xw88XWeM2Qm0Ajbg/3maGdS9FSf2s7o2nQ0ctdPPWvvs4cfGmM/x\nf37zI1fWEfnGmDhr7SFCf06VP8uW+HeGRdL/gDXW2kcrNxzne1DnKv2DM4mg+eOAhvD5AQwBQk7N\nWGtXA6sDj2cZYzKMMZ662OlfOVOMMRH7/mkEH8QYkwKMBS611u6t3GaMmRo0UhoCLI9wfTcaY0YG\nHmcCzfHvUMVauxFINsa0N8Z4gUuBTyJZX6CulkC+tbak0nJjjHnFGOMK1Hc63x0FEmmf8d2vr6uA\nKZXaPwGuBjDG9AG2W2sPRqo4Y8yNQIm19uGq2qv6HkSovrcD+3zA/4955f8P6vXzC9IfWBKqwRjz\na2PM9YHH2UBOHYV7qEyJ2PdPI/jvuxZoCrwRNJ/3Bf7D0N4NjNpnG2MO4d8bHrHRe8Ak4JXAz71o\n4A7gBmNMnrX23cDzVwN9X7fWfhvh+sA/n7j78BNjzAP4j/iYZYzZgn9EVQFMisSOL2NMX+BvQHug\n1BhzNXAj8Lwx5jZgE/BCoO9rwE+stTONMQuMMTMDtd4V4fqaAUXGmK8C3VZaa+88XB8hvgeV/0Gt\n4/qeAl43xhQC+YGa6uXzO0aNw/F/F9dV6vu+tfZy4BVgojHmdvw5eGsdlRcqU24Bxkfi+6frwYuI\nOJSmaEREHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxqP8HDAaTvoqVxrIAAAAASUVO\nRK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f26b949b9b0>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "stream",
          "text": [
            "Binning Statistics: Energy = -3.1408 +- 0.0245, Auto correlation Time = 0.5480\n",
            "Step 20, dE/|E| = 0.1398, elapse = 2.4577\n",
            "Binning Statistics: Energy = -3.1640 +- 0.0257, Auto correlation Time = 0.5925\n",
            "Step 21, dE/|E| = 0.1334, elapse = 2.4249\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "error",
          "ename": "KeyboardInterrupt",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-23-efe984ae8b77>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0mt0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 31\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0menergy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprecision\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlearning_rate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0muse_cuda\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0muse_cuda\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     32\u001b[0m     \u001b[0mt1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Step %d, dE/|E| = %.4f, elapse = %.4f'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0menergy\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mE_exact\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mE_exact\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mt0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-22-28ddb3209194>\u001b[0m in \u001b[0;36mtrain\u001b[0;34m(model, learning_rate, use_cuda)\u001b[0m\n\u001b[1;32m     15\u001b[0m         \u001b[0;31m# as well as the precision of energy.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m         \u001b[0msample_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvmc_sample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_config\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_bath\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m200\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_sample\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m         \u001b[0menergy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0menergy_grad\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprecision\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvmc_measure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlocal_measure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_spin\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m         \u001b[0;31m# update variables using steepest gradient descent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-13-97345af00a71>\u001b[0m in \u001b[0;36mvmc_measure\u001b[0;34m(local_measure, sample_list, measure_step, num_bin)\u001b[0m\n\u001b[1;32m     62\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mmeasure_step\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     63\u001b[0m             \u001b[0;31m# back-propagation is used to get gradients.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 64\u001b[0;31m             \u001b[0menergy_loc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_loc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlocal_measure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     65\u001b[0m             \u001b[0menergy_loc_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0menergy_loc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     66\u001b[0m             \u001b[0mgrad_loc_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrad_loc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-19-80c2aa68884e>\u001b[0m in \u001b[0;36mlocal_measure\u001b[0;34m(self, config)\u001b[0m\n\u001b[1;32m     41\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     42\u001b[0m         \u001b[0;31m# E_{loc}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 43\u001b[0;31m         \u001b[0meloc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menergy_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mansatz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpsi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_numpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpsi_loc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     44\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0meloc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_loc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-21-addfdba35f3c>\u001b[0m in \u001b[0;36mheisenberg_loc\u001b[0;34m(config, psi_func, psi_loc, J)\u001b[0m\n\u001b[1;32m     32\u001b[0m         \u001b[0mconfig_i\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m         \u001b[0mconfig_i\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mflip\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 34\u001b[0;31m         \u001b[0meng_i\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mwi\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mpsi_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig_i\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mpsi_loc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     35\u001b[0m         \u001b[0macc\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0meng_i\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0macc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-19-80c2aa68884e>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m     41\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     42\u001b[0m         \u001b[0;31m# E_{loc}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 43\u001b[0;31m         \u001b[0meloc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menergy_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mansatz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpsi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_numpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpsi_loc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     44\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0meloc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_loc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-20-d39c01e6e263>\u001b[0m in \u001b[0;36mpsi\u001b[0;34m(self, v)\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0mfloat\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mprobability\u001b[0m \u001b[0mof\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m         '''\n\u001b[0;32m---> 33\u001b[0;31m         \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdouble\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     34\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_cuda\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcuda\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mh_bias\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcosh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv_bias\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
          ]
        }
      ]
    },
    {
      "metadata": {
        "id": "t1F2-L3F1rBr",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Challenge\n",
        "J1-J2 model is an advanced version of Heisenberg model that takes 2nd nearest neighbor coupling terms into considerations. Its hamiltonian is \n",
        "\n",
        "$$H=J_1\\sum\\limits_i \\vec S_i\\cdot \\vec S_{i+1}+ J_2\\sum_i\\vec S_i\\cdot \\vec S_{i+2}$$\n",
        "\n",
        "The sign of its ground state wave functions $\\psi(x)=\\langle x|\\psi\\rangle$ can be either positive or negative. Since real valued RBM is not able to describe such a state.\n",
        "\n",
        "With the following local energy function, write your own ansatz to obtain its ground state.\n",
        "\n",
        "As a start point, you can try with 20 sites, $J_2=0.2$.\n",
        "Obtain the best precision as you can.\n",
        "\n",
        "Below is a table of DMRG results for benchmark, which is believe to be \"exact\" for 1D systems.\n",
        "\n",
        "| &nbsp;&nbsp;&nbsp;$J_2$ &nbsp;&nbsp;&nbsp;  | 20             | 30             | 40             | 100            |\n",
        "| -------------------- | -------------- | -------------- | -------------- | -------------- |\n",
        "| 0.                   | -8.90438652988&nbsp;&nbsp;&nbsp; | -13.3219630586&nbsp;&nbsp;&nbsp; | -17.7465227719&nbsp;&nbsp;&nbsp; | -44.3229467082&nbsp;&nbsp;&nbsp; |\n",
        "| 0.2                  | -8.20291625218 | -12.2770471706 | -16.3566615295 | -40.8572924302 |\n",
        "| 0.5                  | -7.5           | -11.25         | -15            | 37.5           |\n",
        "| 0.8                  | -8.46127240196 | -12.6588455544 | -16.8706800952 | -42.0700632095 |\n"
      ]
    },
    {
      "metadata": {
        "id": "QkR3KNCw1wFx",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        "def J1J2_loc(config, psi_func, psi_loc, J1, J1z, J2, J2z, periodic):\n",
        "    '''\n",
        "    local energy for 1D J1, J2 chain.\n",
        "\n",
        "    Args:\n",
        "        config (1darray): bit string as spin configuration.\n",
        "        psi_func (func): wave function.\n",
        "        psi_loc (number): wave function projected on configuration <config|psi>.\n",
        "        J1, J1z, J2, J2z (float): coupling strenghs for nearest neightbor, nearest neighbor z-direction, 2nd nearest neighbor and 2nd nearest neighbor in z-direction.\n",
        "        periodic (bool): boundary condition.\n",
        "\n",
        "    Returns:\n",
        "        number: local energy.\n",
        "    '''\n",
        "    nsite = len(config)\n",
        "    wl, flips = [], []\n",
        "    Js = [J1, J2]\n",
        "    Jzs = [J1z, J2z]\n",
        "    for INB, (J,Jz) in enumerate(zip(Js, Jzs)):\n",
        "        # J1(SzSz) terms.\n",
        "        nn_par = np.roll(config, -INB-1) * config\n",
        "        if not periodic:\n",
        "            nn_par = nn_par[:-INB-1]\n",
        "        wl.append(Jz / 4. * (nn_par).sum(axis=-1))\n",
        "        flips.append(np.array([], dtype='int64'))\n",
        "\n",
        "        # J1(SxSx) and J1(SySy) terms.\n",
        "        mask = nn_par != 1\n",
        "        i = np.where(mask)[0]\n",
        "        j = (i + INB+1) % nsite\n",
        "\n",
        "        wl += [J / 2.] * len(i)\n",
        "        flips.extend(zip(i, j))\n",
        "        \n",
        "    # calculate local energy <\\psi|H|x>/<\\psi|x>\n",
        "    acc = 0\n",
        "    for wi, flip in zip(wl, flips):\n",
        "        config_i = config.copy()\n",
        "        config_i[list(flip)] *= -1\n",
        "        eng_i = wi * psi_func(config_i) / psi_loc\n",
        "        acc += eng_i\n",
        "    return acc"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "qRBaeQgTCgb9",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## Hints\n",
        "* spliting phase and amplitude parts of a wave function is suggested in\n",
        "\n",
        "    `Cai, Zi, and Jinguo Liu. \"Approximating quantum many-body wave functions using artificial neural networks.\" Physical Review B 97.3 (2018): 035116.`\n",
        "* using **symmetries** like invariance under translation ,spin flip and inversion in your model to make training more stable.\n",
        "Ground states of J1-J2 models with different $J_2$ or $N$ lies in different symmetry blocks.\n",
        "* use stochastic reconfiguration or Adam as a better **optimization** method.\n",
        "* rewrite the part of code to get local energy using C/C++/Fortran for speed up.\n",
        "* use tensorflow or your own modules to support **complex values**.\n",
        "* for large scale 2D J1-J2 model, its ground state is an **open question yet**!"
      ]
    },
    {
      "metadata": {
        "id": "SKjnfSMq6AWf",
        "colab_type": "code",
        "colab": {
          "autoexec": {
            "startup": false,
            "wait_interval": 0
          }
        }
      },
      "cell_type": "code",
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}